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Abstract: The next-to-lcading order (NLO) evolution of the parton dis- 
tribution functions (PDF's) in QCD is the "industry standard" in the 
lepton-hadron and hadron-hadron collider data analysis. The standard 
NLO DGLAP evolution is formulated for inclusive (integrated) PDFs and 
is done using inclusive NLO kernels. We report here on the ongoing project, 
called KRKMC, in which NLO DGLAP evolution is performed for the ex- 
clusive multiparton (fully unintegrated) distributions (ePDF's) with the 
help of the exclusive kernels. These kernels are calculated within the two- 
parton phase space for bremsstrahlung subset of the Feynman diagrams of 
the non-singlet evolution, using Curci-Furmanski-Petronzio factorization 
scheme. The multiparton distribution with multiple use of the exclusive 
NLO kernels is implemented in the Monte Carlo program simulating multi- 
gluon emission from single quark emitter. With high statistics tests (~ 10 9 
events) it is shown that the new scheme works perfectly well in practice 
and is equivalent at the inclusive level with the traditional inclusive NLO 
DGLAP evolution. Once completed, this Monte Carlo module is aimed as 
a building block for the NLO parton shower Monte Carlo, for W/Z pro- 
duction at LHC and for ep scattering, as well as a starting point for other 
perturbative QCD based Monte Carlo projects. 
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1. Introduction 

The so called next-to-leading-order (NLO) parton shower Monte Carlo 
(MC), usually referred to as a highly desirable type of the calculation tool for 
perturbative QCD predictions at LHC is commonly believed to be unfeasible 
in practice. The ongoing project presented here demonstrates proof of the 
existence of such a NLO parton shower MC for the initial state QCD, albeit 
for a limited subset of diagrams of the non-singlet NLO QCD evolution of the 
parton distributions functions (PDFs) The project involves re-calculation 
of the NLO evolution kernels in the exclusive (fully unintegrated) form, 
following Curci-Furmanski-Petronzio method |T] for the DGLAP type [2] 
of the QCD evolution of PDFs. A prototype MC with the new exclusive 
kernels performs exactly the NLO DGLAP evolution of PDFs on its own 
(no external pretabulated PDFs needed!), following closely the collinear 
factorization theorems of QCD and the standard "Matrix Element x Phase 
Space" approach. Once completed, this project will open new avenues for 
a new class of perturbative QCD calculations in form of Monte Carlo event 
generators for the coming two decades of the LHC experiments. 

What are the main aims and assumptions of the project? The new NLO 
Parton Shower Monte Carlo for QCD initial state radiation should be: 

• based firmly on Feynman diagrams (matrix element) defined within 
the standard Lorentz invariant phase space, 

• based rigorously on the collinear factorization of the classic works, 
see or [I], 

• implementing exactly NLO DGLAP evolution in MS scheme, 

• implementing the exclusive PDF (ePDF) defined within the standard 
multiparton Lorentz invariant phase space^ 

• performing NLO evolution by the MC itself, with help of the new 
exclusive NLO kernels. 

The above priorities can be compared with ref. [6], where construction of 
NLO parton shower MC is also advocated. 

In the realization of the presented KRKMC project very strong emphasis 
is put from the very beginning on testing all theoretical and practical ideas 
with series of high quality numerical tests. In fact, every major milestone 
of the project realization is marked by the construction and testing the 
corresponding numerical MC prototype. In this contribution we shall report 
results form first few MC exercises of this kind. Unfortunately, we shall 



1 Similar kind of PDF is also referred to as a "fully unintegrated" PDF [5] . 
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not be able to give all the details of these MC constructions in this short 
presentation - we hope to report them in detail in separate publications 
soon. 

Let us attempt to describe very briefly the basic features of the theoret- 
ical QCD model of the parton distributions in the KRKMC NLO MC. In 
the folowing formula implemented in the KRKMC program, 

°° f n d 3 ki 

D{x,Q) = ^2 Tl^W P (n) (P;h,k 2 ...,k n ) & S (k u k 2 ...,k n )<Q S 1-x=^ch, 
n =o J i=i LK i 

< , (1) 

the exclusive parton distribution function, ePDF, is the integrand p^ '. In 
the above we denote the lightcone variables of the emitted partons as a% = 

W = ki 2E i ' wri ere fcf = (k®, kj, kf ), i = l,..,n are 4-momenta of these 
partons, while the initial parton 4-momentum is qo = P = {E, 0,0, E). For 
the variable S (factorization scale) we use one of the following variables 



Sq 



i=l 



1/2 ^kil Ikol Ik 



S r = max 



1| |K 2 | 



|k j (2) 

S v = max ( -^=,, . . . ) , S T = max(|ki|, |k 2 |, . . . , |k n |), 
la\ Ja 2 Ja n J 



S q = Q is used in ref. [T], while in the following we shall usually work 
with maximum transverse momentum St- Variable S r related to maximum 
rapidity of the emitted real partons is generally our preferred choice and 
occasionally we shall use maximum lightcone minus variable S v . 

The integrands in eq. ([I]) are finite due to built-in regulators of the 
collinear divergences. The actual integrands p^ in eq. |l| are the result of 
projecting to the LO+NLO level the differential distribution 

p^(P; k 1 ,k 2 .-., k n ) -» p^(P; h,k 2 ...,k n ), (3) 

where p^ is coming from the UV subtracted Feynman diagrams^] using 
methodology of the "factorization theorems" of QCD in the physical gauge, 
such that the inclusive D(x, Q) obeys exactly the standard evolution equa- 
tion 

j^D(x, Q) = j dzdu y LO+NLO (Q, z) D(u, Q) 5 X=ZU (4) 



of the NLO DGLAP with the standard inclusive kernels of the MS scheme 

2 

^LO+NLO i 



3>^+^(Q, z) = °^-p(°\z) + ( PM(Z). (5) 

Z7T \ Z7T J 



It includes IR regulators of the collinear and soft singularities. 
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How to perform the projection of eq. Q is of course the main theoretical 
issue in the project. In particular even the existence of this projection, such 
that eq. Q is fulfilled, is highly nontrivial and still open question. 

On one hand, we insist very strongly on the perfect compatibility of our 
ePDF's with the NLO DGLAP in MS scheme in the sense of eq. Q and we 
believe that it is feasible. On the other hand, this requirement should not be 
overstressed, as we do foresee in the next step going beyond NLO DGLAP 
in the exclusive form, for instance towards BFKL. However, in our opinion 
such an extension has to wait until the exclusive NLO DGLAP evolution in 
the MC form is well established and tested. 

The multiparton distribution p^ describes the chain of the parton emis- 
sions out of the single emitter parton coming from the hadron beam and 
absorbed in the hard process. In the classic LO parton shower MC p( n ) is 
rather simple, just the product of the LO kernels times Sudakov formfactor. 
Momenta of the MC event are generated according to LO level p^ n \ typ- 
ically using backward evolution algorithm of ref. [7] , or the new constrained 
MC algorithms of refs. [H |9j. On the other hand, the new and more com- 
plicated p( n ) at the NLO level will be basically the convolution of the 2PI 
kernels of the collinear factorization scheme of refs. [3] [JJ, truncated at the 
complete NLO level (in the axial gauge). Of course, p^ has to be defined 
in four dimensions, d = 4, and the removal the IR dimensional regulator 
from the real emission part of the Curci-Furmanski-Petronzio scheme [JJ, 
without spoiling main features of the standard MS scheme, is one of the 
main practical issues. 

Since the main goal of the KRKMC project is to reproduce NLO DGLAP 
evolution in QCD in the exclusive way, in the Monte Carlo, while its stan- 
dard inclusive version is a well established technique and serves well to 
describe lepton-hadron and hadron-hadron collisions, it is therefore a valid 
question to ask: "why bother to do the same in more complicated way". 
The answer is that it is worth to redo it in the exclusive MC form, because 
there are numerous potential gains once it is available; it may serve either 
as a building block in a bigger new MC project, or as an excellent starting 
point for many new types of the MC parton shower type MC projects. Let 
us list some of the possible applications: 

• Extension towards CCFM [ID], BFKL [IT] (low x limit in PDFs), with 
the correct soft limit and built-in colour coherence. 

• More realistic description of the heavy quark thresholds in PDFs and 
the hard process. 

• The use of exact amplitudes for multiparton emission in the soft and 
collinear limits - the analog of Coherent Exclusive Exponentiation in 
QEDp]2[l3]. 
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• Better connection between hard process ME and the shower parts, as 
compared with MC@NLO jH] and the other similar approaches of this 
class. In particular we expect no negative weight events, no ambiguity 
of defining last emission before hard process, etc. 

• Providing better tool for exploiting HERA DATA for LHC, in partic- 
ular for fitting i*2 directly with the MC programs. 




Fig. 1. Feynman diagrams contributing to the calculation of the non-singlet NLO 
kernel, except two-loop self-energy graphs. 

2. Constructing and testing Exclusive kernels 

The first inevitable step in the construction of the multiparton MC im- 
plementing complete NLO DGLAP in the exclusive way is to construct and 
test the exclusive NLO evolution kernels within the two parton phase space 
in d = 4, which later on will be used many times in the n-parton emission 
ladder. Since our ambition is to be fully compatible with the NLO DGLAP, 
it is therefore natural that we take as a starting point the calculation of 
the NLO kernels in the diagrammatic way according to Curci-Furmanski- 
Petronzio (CFP) scheme [I], rather than technique of calculating anomalous 
dimensions of the Wilson operators |15|, I16j . In the CFP scheme one ex- 
ploits the statement of ref. [3] (EGMPR) that in the physical (axial) gauge 
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all collinear singularities are located on the rungs of the ladder diagram 
between the 2PI kernels. Feynman diagrams contributing to non-singlet 
NLO kernels are depicted in Fig. [T] In the following discussion we will limit 
ourselves to two bremsstrahlung diagrams inside the dashed box marked in 
this figure, in fact to C\ part of them. 

2.1. Re- Calculating NLO DGLAP kernels following CFP 

Recalculation of the NLO DGLAP kernels is mandatory because we 
need access to the two parton differential distributions in the internal phase 
space of the kernels. Moreover, we shall need explicit expressions for the vir- 
tual contributions to NLO DGLAP kernels from virtual Feynman diagrams, 
which are not available in ref . [I] , nor in the papers |17[ [18] , where the cal- 
culations of refs. [TJ [19] were reproduced and cross-checked once again. 

Collinear factorization theorem of EGMPR |3j improved and customized 
to the use of MS dimensional regularization scheme by CFP [1] states that 



1 f Q 2 \ ( 1 



G 







l-(l-P).Jjf j } 1 _( PKo .^i 



1-(1-P)-K )) (6) 



r(a,-) = ( — — — I = 1 + K + K&K + K® K®K + 



ej \1-K 

K = PK0 'l-(l-P)-K > C = C °-1-(1-P).K - 

where Kq is 2-particle irreducible (2PI) kernel, which is free of the collinear 
divergences. The ladder part T corresponds to the MC parton shower and C 
is the hard process part. The projection operator of CFP scheme is defined 
as follows 

P — Pspin Pkin PP (7) 

and it consists of the kinematic projection operator Pfcj n , spin projection 
(averaging) operator P sp i n and the pole part PP extracting ^ , k > partrl 
Multiplication symbol • denotes full phase space integration d n k (the only 
source of collinear divergences), while convolution <S> involves the integration 
over the 1-dimensional lightcone variable x only. 



3 The above formula should be treated with special care, see ref. [J], due to non- 
associativity of the products of the kernels with insertions of the (1 — P) operator. 
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In the CFP methodology the NLO kernel is extracted from the second 
order term in the expansion: 

r = — ?— = 1 + PK + PK • [(1 - P) • K ] + [PK ] ® [PK ] + • • • , 
1 — A 

C = C °1-(1-P)K 

= Co + Co ■ (1 - P) • K + Co • (1 - P) • K ■ [(1 -P)-K )} + ... 

(8) 

Following the "pole-part method" of CFP the NLO kernel is obtained as 
followfS 

= ^P(a,x) = ^P {0) (x) + (^) 2 pW(x) = a ^-Res x T(a,x) 
2ir 2tt \2ir/ oa 

= Re Sl {PK } + 2Re Sl {PK ■ [(1 - P) - K }]. 




Fig. 2. Extraction of NLO kernel for bremsstrahlung diagrams. 

The role of the P projection operator in the CFP scheme in the case 
of our bremsstrahlung diagrams is illustrated in Fig. |2j The upper line in 
Fig. [2] shows the three Feynman diagrams, defined in the standard phase 
space which is closed from above by the operator P. It employs also cut-off 
parameter Q = [ip to limit phase space from the above. The well known 



4 Resi denotes coefficient in front of - pole in the Laurent expansion around e = 0. 
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property of the physical gauge [20] tells us that 2nd diagram in Fig. [2] 
contributes double collinear log (or 1/e 2 ), while the 1st and 3rd only single 
collinear log (or The insertion of P+(l— P) in between two gluon rungs 
in the 2nd diagram changes nothing in the total sum, however, the part 
containing (1 — P) is now only single logarithmic, and in the factorization 
procedure it is combined with the last "crossed diagram" , forming together 
the NLO contribution to evolution kernel, as depicted in the lower part 
of Fig. [2} Let us remark that the remaining part of 2nd diagram with two 
gluon rungs and two P's is just the square of the 1st diagram, so it represents 
the "idealized" 2nd order LO contribution (the beginning of the LO series). 
The NLO part (dashed box) is integrated over the 2-gluon phase spacej^jand 
associated with the single gluon 1st diagram (see dashed arrow), forming 
the NLO part of the standard inclusive DGLAP kernel. We want to stress 
that this procedure combines 1 and 2 real gluon contributions (plus virtual 
ones) into a single inclusive object, the standard NLO inclusive kernel. 

Our main aim is to recover the internal 2-gluon differential distribution 
of the NLO kernel in the MC event generator. In a sense our aim is to undo 
the procedure marked by the dashed arrow in Fig. [I| Let us now concentrate 
on the above internal 2-gluon differential distribution of the NLO kernel for 
the bremsstrahlung diagrams (Cp part). 



2.2. Two-gluon differential distribution inside NLO DGLAP kernel 




Fig. 3. Kinematics of the 2 gluon emission. 



Before we describe the 2-gluon differential distribution let us define kine- 



5 Constraining total loss of the lightcone variable. 
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Xi 



C • q% 
C • <?o' 

o, p 



C 7 ^ ~ 2E' 



q = (E,0,0,E), 



C = (1,0,0,-1), c 2 = o, 

kf = 0, k{ = (k®, kj, kf ). 



qo 



(10) 



The second order contribution to T, for the 2nd bremsstrahlung diagram 
in Fig. [2] in d = 4 + 2e dimensions, reads 



Co^o^o = C 2 a 2 / d0 ^ da ^8 1 ^ x=ai+a2 I d 2+2e k 1 d 2+2£ k 2 /x" 4e e(5(fci,fc 2 ) < Q) 



j ' Ti(qi,Q!2) , r 2 (ai,a 2 ,e) k| , r 3 (ai,a 2 ) 2kj ■ k 2 



aia 2 



a 2 



(ii) 



with the notation explained below and the extracted kernel being indepen- 
dent of the actual choice of Co in eq. (|ll]Q The C 2 F part of the NLO kernel 
extracted from V function for all 3 bremsstrahlung diagrams including 2 soft 
counterterms reads 



^ N ( x ) = 2{ J ~2k~^ J ~2k^ < ^ 1=max (l k il'l k 2l)/ ( 9 fii-x=ai+a2 b%(ki,k2), (12) 



We use Co = Pspin0s{ki,k 2 ,~-,kn)<Q> wnere S is one of these in eqs. (|2| 
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where 

bg(k U k2 



tLadd. 



(ki,k 2 ) 



+ T 2 v (a 1 ,a 2 

b Count \k x ,h 
- q 2 {h,k 2 ) 



(aC F ) 2 


^Ladd 


16(2vr) 2 




1 


fe) { 


q A (h 
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^2){ 


q 4 (ki 





■b 0ount (k 1 ,k 2 ) 

b Count ik 2 ,k 1 ) + b XLad -(k 1 ,k 2 ) 



Ti(ai,a 2 ) r 2 (ai,a 2 ,0) kj T 3 (ai,a 2 ) 2ki • k 2 
a±a 2 

2T?(ai,a 2 ) 



On 



C>2 



ai« 2 

2ki-k 2 , 2ki-k 2 , rrXf .(kx-k,) 2 
i -3- + r 2 (a 2 ,ai) ry 



k 2 , T 2 (ai,a 2 ,0) 1 



+ 2f(ai,a 2 )- 
T 2 (ai,a 2 ,0) 



k 2 k 2 

K 1 K 2 



9 4 (0,A: 2 



k 2^<k^ (l_ ai )2 k 2 k 2-kf<k^ 



k? + -k^ + 2ki k 2 . 



Q'2 



(13) 



where 7-trace factors and the second quark propagator are 

Ti(a\,a 2 ) = (1 + x 2 + x 2 )aia 2 , T 3 (ai, a 2 ) = (1 + x 2 + x 2 )xi, 
T 2 (ai,a 2 ,e) = (1 + x\){x 2 + x\) + eT 2 (ai,a 2 ), 
I^(ai, a 2 ) = (1 - x x ) 2 {x 2 + x\) + (x- xi) 2 (l + xj), 
T£(a\, a 2 ) = x(l + x 2 - «ia 2 ), If (ai, a 2 ) = 2(1 + x 2 ), 
T|(ai, a 2 ) = a?(l - ai) + (1 + x 2 )(l - a 2 ), 



(14) 



The integral (13) is already in d = 4 dimensions and it is finite. For sim- 
plicity we have omitted ternj^] T' 2 from the e part in the 7-trace leading to 
T 2 . 



The integrations in eqs. (12) have been performed analytically: 

2 



y N ( z ) 



r 2 a 



~ / 1 + 3x 2 , 2 



irJ \16{l-x) 



x 3 
ln^(a;) H ; — ln(x) H — (1 — x] 



(15) 



reproducing the result of CFP pp. On the other hand the distribution 
b 2 (k\, k 2 ) of eq. (13) has been plugged into the general purpose MC inte- 
grator/simulator FOAM and the result of numerical MC integration was 
compared with formula of eq. (15), see Figj4j Both NLO results agree to 

7 We keep them in the actual program and in the following numerical tests. 
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O.S 0.9 




Fig. 4. Comparison of numerical MC and analytical integration of (fci, k 2 ) con- 
tribution to exclusive NLO kernel. 



within statistical MC relative error of a couple percent. The same error 
with respect to LO+NLO is only ~ 1CT 3 . 

However, our real interest is not the integral but rather the integrand 
itself, as it will be used as a building block in the multiparton distribution 
for the MC, hence in the following we shall analyze the shape of b^(k\, fa) 
and its singularity structure in a detail. 



2.3. Analyzing internal structure of the NLO DGLAP kernel 

The activity of the previous section, culminated in Figj4] was just the 
warm-up exercise, in which we have checked whether we are able to build a 
numerical model for the NLO kernel in the exclusive form. The next step, 
which we pursue in the following, is to examine the shape and singularity 
structure of the integrand. The final step will be to implement the distri- 
bution (k\, fa) several times in the actual multiparton distribution of the 
NLO parton shower MC program, see next Section. The extension of the 
present analysis of the IR singularity structure to more Feynman diagrams 
and genuine non-abelian contributions ~ CfCa is reported in ref. [21] in 
these proceedings. 

To begin with we check how does the integrand (fa, fa) look like in the 
Sudakov-type variables ln(A;f /k^) and ln(c>!i/a:2), keeping two constraints: 
max(fc^ / fcj) = Q and a\ + 02 = 1 — x. In particular we are interested in 
the soft limit of one of the two gluons, for instance k\ —> and 02 — ► 0. 
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Fig. 5. Contribution to unintegrated NLO kernel from one of the double- 
bremsstrahlung diagrams (top left-hand plot). Subtraction of the internal sin- 
gularity, k% — > 0, is done using soft counterterm (top right-hand side). Subtracted 
distribution is also shown (bottom). Lightcone variable is fixed, x — 0.3. 



The single diagram of bremsstrahlung with the soft subtraction 

b Ladd ik2,h) -6 Coun *-(fc 2 ,A;i) 

is visualized in Fig. [5] In this figure the upper-left-hand plot represents the 
distribution from single Feynman diagram of the double bremsstrahlung 
type. The "cliff" along the line of equal virtualities of the emitted gluons 
ki/y/ax = /c<f/y / 02 represents well know "ordering effect" in the physical 
gauge due to propagator of the emitter quark. The soft counterterm of 
CFP depicted in the upper-right-hand plot removes soft singularity in the 
limit k\ — > 0. However, the remaining structure near kf ~ k 2 features 
very strong cancellations between two triangular regions of the doubly- 
logarithmic size. Such structures, if they were really present, would render 
any MC implementation of the NLO unintegrated kernel unfeasible due to 
strong MC weight fluctuations. 

Fortunately, the inclusion of the second bremsstrahlung diagram 

b Ladd -(ki,k 2 ) - b Count -{k u k 2 ) + b Ladd -(k 2 , ki) - b Count (k 2 , k x ) 
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Fig. 6. Contribution to unintegrated NLO kernel from two double-bremsstrahlung 
diagrams with subtraction, x = 0.3. 




Fig. 7. Contribution to unintegrated NLO kernel from all three double- 
bremsstrahlung diagrams with subtraction, x = 0.3. 

makes the distribution more friendly for the future MC use, as shown in 
Fig. [6j Nevertheless the single-log structure along the line j \fot\ = 
&2 '/ \foL2 seems to contradict the eikonalization of the bremsstrahlung dis- 
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tributions in the soft limit [22] and would cause the problems for the MC 
use. The size of the remaining uncanceled part is much smaller than the LO 
distributions. In this plot we have used x = 1 — ct\ — oli = 0.3 and the inte- 
gration over azimuthal angles of the gluons has been performed. However, 
we have checked that general features of the above results remain the same 
for any values of x and any gluon azimuthal angle. In ref. [21] similar plot 
is shown for b Ladd "(ki,k2) + b Ladd \k2,k\), that is for the bremsstrahlung 
matrix element without soft counterterm subtraction. 

In Fig. [7] the above single-logarithmic structure disappears when finally 
the "crossed" diagram (interference) is included. As we see, the interference 
diagram (see upper-right-hand plot) features the same single-log structure 
along the line of equal virtualities, which corrects the soft eikonal limit and 
the sum is nonzero only in the region of kf ~ k% and ai ~ ^{1 — x) ~ 1. In 
other words, (ki, &2) plays a role of the "short range correlation" function 
on the Sudakov plane which enters into the game only when both gluons are 
not soft and they are not in the so-called "strong ordering" regime. Such 
a behaviour looks also very friendly for the Monte Carlo use, if we plan to 
use b^ (k\, /C2) as a component in the MC weight correcting LO distribution 
to the NLO level. 

3. Re-insertion of exclusive NLO kernels into LO Monte Carlo 

Once exclusive representation of the NLO evolution kernel in two parton 
phase space is constructed and analyzed, the next question is how to use it to 
construct the multiparton distribution of the PS MC. We shall show in the 
following how to solve this problem for the Cj, part of the bremsstrahlung 
diagrams contributing to the non-singlet NLO kernel. 

3.1. General scheme of exclusive NLO insertion 

The very essence of the general methodology of the multiple "insertion" 
of the exclusive NLO kernel in the ladder of CFP scheme is depicted in 
Fig. [8} for the NLO MC distributions with up to four gluons. The extension 
to more gluons is straightforward. In Fig. [8] terms with boxes marked as 
LO only, form the traditional LO distributions being the product of the LO 
kernels, times Sudakov formfactor resuming virtual corrections, see below 
for more details. Note that LO box includes trivial single real parton phase 
space. On the other hand, the NLO box, which shows up for the first time 
for 2 real gluons and is defined in the top partj^] of Fig. |8j features two real 
parton phase space. The case of 2 gluons in Fig. [8] is trivial, because it is 

8 Of course, we keep in mind that NLO kernels include virtual contribution within the 
single real parton phase space, similarly as LO box. 



IFJPAN-IV-09-3 printed on May 9, 2009 



15 



equivalent to the definition of the NLO kernel. The first nontrivial term in 
the decomposition of Fig. [8] is that for 3 real gluons, where the two-gluon 
NLO part inserted either before or after the LO kernel. In case of four gluons 
in Fig. [8] we have 3 terms with one NLO insertion in 3 possible ways and a 
double NLO insertion for the first time, albeit without any LO spectators. 
The extention of this pattern to any number of real gluons is not difficult. 



— — ~ — — — — — — — ~ — — 

Fig. 8. General scheme of exclusive NLO insertions. 




} 



The real problem is how to define, or rather deduce from the original 
matrix element, the convolution of many LO and NLO factors within the 
standard Lorentz invariant phase space, such that the infinite sum of such 
convolutions obeys exactly the LO+NLO evolution equation of eq. Q. 

At first glance the above task may seem unfeasible. Nevertheless, we 
shall show in the following how to do it, and demonstrate the first Monte 
Carlo implementation, albeit only for the bremsstrahlung subset of the di- 
agrams of the previous Sections. Not surprisingly, the solution is found by 
examining carefully the original ra-parton matrix element coming directly 
from the Feynman diagrams, and the procedures of the QCD factorization 
theorems used to extract the "idealized" LO+NLO part from it. 



3.2. NLO insertion for DGLAP evolution with inclusive kernels 

We start the above difficult task with the warm-up exercise in which 
the NLO insertions are defined within the traditional iterative solution of 
the DGLAP evolution equation with LO+NLO inclusive kernels. The cor- 
responding simple Markovian MC program with the inclusive LO+NLO 
kernels will also serve as a valuable testing tool of a more complicated new 
parton shower MC with the exclusive LO and NLO kernels. 

Changing slightly notation, the traditional NLO DGLAP evolution reads 
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d t D(t,x) = a D(t))(x) = a [ dzdx 5 x=zy l J'{z)D(t,y), (16) 

Jo 

where the inclusive LO+NLO kernel is 9(z) = 9^°\z) + a?^\z), a = C F % , 
t = In Q = In /i and a is non-running. The NLO part of the kernel 
a 2t ?( l \z) = y N (z) is exactly that of eq. (12). The iterative solutions of 
the above equation reads as follows^ 



00 / n ft \ 
D(t,x) = e- s 5 x=1 + e~ s Y, II / dt * (^e\ z i) + « 2 ^ (1) W) <Un*> 

n=l \i=lJ t *-l ) 

^ \z)=(-^- s + l)s x = 1 + yf\z), ?£\z)- 1 + z2 



2(1 -z) 



n-z>6, 



S-a(l- lt)) [hi~-~]. 



(17) 



In the Markovian MC exercise parton multiplicity n and the entire MC 
event {U,Zi,i = 1,2,..., re} is generated, using only LO kernel 9<®(z). The 
NLO content is introduced using MC weight 

nr=i(^%) + _ * T{d) rgu gg(g) (18) 



9 Virtual LO corrections regularizing soft divergence 1/(1 — z) are resummed into Su- 
dakov formfactor exp(— S). The NLO part in S is set to zero and J 7(z)dz 7^ 0. 
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where {d} is the set of all 2 n partitions d = (d±, cfo, d^, d n ), di = 0,1. 
r (d) = Yli di is the number of NLO factors in a given partition. The 
decomposition in Fig. [8] is essentially that of the numerator in the above 
expression truncated to terms of the order a n , n < 4. 

Fig. [9] presents result from the MC run for evolution from Qo = e*° = 
lGeV to Q = e* = lOOOGeV and a = 0.2. In the plot we see LO result for 
r = 0, the NLO result, and the slices with the defined number of the NLO 
insertions being r = 0,1,2,3,4. 

The important conclusion from this introductory exercise is that for 
any practical purpose we may truncate the surrj^l over NLO insertions to 
r(d) < 3 or even r{d) < 2. 



3.3. Introductory step: extraction of LO from the exact matrix element 




Fig. 10. Extraction of LO from the exact matrix element for 3 real gluons. 



As already stressed, the main problem in devising LO+NLO exclusive 
distribution is how to combine NLO insertions and LO parts. We have 
solved it by examining once again how LO+NLO distributions are extracted 
from the exact multiparton matrix element (i.e. from Feynman diagrams). 



For the subset of the bremsstrahlung diagrams under discussion. 



18 



IFJPAN-IV-09-3 printed on May 9, 2009 



This cannot be done without careful re-examination of this procedure for n 
partons at the seemingly trivial LO level. 

Let us do it for n = 3, that is for the emission of 3 gluons out of an 



emitter quark. The whole procedure is depicted in Fig. 10 and it goes in 
two steps: we start with re! = 6 diagrams squared defined in the entire 
phase space (and 1/nl Bose factor in front). The explicit squared sum over 
re! permutations {vr} = (123), (213), (132), (231), (312), (321) is expanded 
into (n!) 2 = 36 terms. 

In the axial gauge in the LO approximation it is allowed to drop all 
interference terms [20], so in the 1st step we do it retaining out of (re!) 2 terms 
only n! = 6 diagonal terms depicted in the second line in Fig. [TUJ Each of the 
n! diagrams covers approximately one simplex v > v ns > v n2 > v ni > vq, 
where Vi = kj / ' ^foi = k~ . This for example is illustrated in the upper 
left plot in Fig. [5] However, as clearly seen in this figure, the shape of 
the distribution is fuzzy at the border of the simplices v\ ~ V2- In other 
words the bremsstrahlung matrix element (squared) of the single Feynman 
diagram in the axial gauge features what we call fuzzy ordering. 

In the 2nd step the insertion of the projection operator P replaces the 
fuzzy ordering (at the simplice boundaries) by the sharp one. In fact the 
MS scheme dictates the LO distribution to be sharp along the lines of equal 
transverse momentum^ Q > fej 3 > k^ 2 > k^ > Qo, as also seen in Fig. [5] 
for n = 2. 

The important lesson from the above LO considerations is that it is 
wise to keep the complete phase space, without any artificial division into 
simplices suggested by the so called ordering property of LO-approximated 
matrix element in the variables like Vi or kf, but instead to associate such a 
division only with the properties of the particular diagram. At the LO level 
there is, of course, one-to-one correspondence among n! squared diagrams 
and n! simplices in k T in the phase space. This may easily be the source of 
the confusion, because this one-to-one correspondence exists only at the LO 
level and breaks down beyond the LO approximation, when one gradually 
goes back towards exact matrix element, for instance by recovering terms 
contributing at the NLO level, or replacing sharp ordering of LO by the 
fuzzy ordering of the exact matrix element. This we have already seen 
at work when studying properties of the NLO unintegrated kernel in the 
previous section. 



11 In the soft limit cii — > the exact matrix element, summed over n! diagrams, is not 
able to tell us whether the LO distribution (also summed up over n\ diagrams) with 
sharp ordering using kj or vt variables is better, as they all coincide in this limit. 
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Fig. 11. Additional NLO contributions forming single NLO insertion (dashed el- 
lipsis) for n — 3 gluon emissions. 



3.4- NLO insertion and LO spectator for n = 3 



The necessity of combining LO kernel with NLO insertion is encountered 
for the first time for 3 gluons. Let us analyze the problem and find solution 
in this simplest beyond-the-LO case. 

Following similar factorization procedure for 3 gluons as in the LO dis- 
cussion of the previous section, let us now retain the contributions of the 
NLO class in addition to the product of LO kernels. The new NLO objects 
are of two kinds: (i) the interference contributions among two adjacent 
emissions in the Feynman diagram, which were neglected in step 1 of the 



previous LO example, see (al) and (bl) in Fig. 11, and (ii) the contribu- 
tion due to restoring of the "fuzzy ordering" of the exact matrix element 
(removing P operator responsible for the sharp ordering) for the same two 



adjacent emissions, see (a2) and (b2) in Fig. 11 



The extra terms of two kinds show up for 3 gluons in two possible loca- 



tions in the emission ladder, as shown in Fig. 11 and form in a natural way 



a single NLO insertion in the presence of a single LO spectator, (al)+(a2) 



and (bl)+(b2) in Fig. 11 As seen in Fig. 11 the LO spectator can be posi- 
tioned either before or after NLO insertion. Let us now translate the above 
description into formulas. 

First of all, we upgrade slightly our LO MC model, by means of adding 
into the game the azimuthal angle of the emitted partons and restoring four 
momenta of all partons. This is done with the help of the usual mapping 
of the evolution variables into four momenta: kj = e**, k~j~ = 2(xj-i — 

Xj)Eh, kj = (kj) 2 /kj~. The distribution of the LO MC for n = 3 emissions 



takes now the following shape within the standard phase space, cf. eq. (17), 
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n L( \ _ -s f d 3 h f d 3 k 2 f d 3 kx -fj 
D 3 (t,X)-e J J J -y-Q di- x=ai + a2+a3 J_|f e t>[k i |>e*o 



p L {h\x. 



P (h\x 2 ) P (AfefcOp (fci|x ) 6*|k 3 |>|k 2 |>|ki| 5 
a 1 + zf 



i-i/ 



16vr (1 - ^)|ki| 2 " 

(19) 

Remembering lessons from the discussion of the LO factorization we rewrite 
the above with the explicit sum over Feynman diagrams squared^] 

r \ _ -S [ ^3 f ^ 3 ^2 f -A- 

D 3 [t,x) — e J —q J —q J -^p- 0i- x=ai+a2+a3 J_J_^>| k ,|>e t o 

3 2 1 i=l 

o[ P L ( k ^\ x 2) P L ( k ^ 2 \ x i) p L (fa 1 \x Q ) 6 , |k, 3 |>|k, 2 |>|k, 1 |, 

' {*} 

(20) 

where fcf = (fc-\ kj, fc 3 ), cti = (&° + kf)/(2E h ), x\ = xo - a^^x^ = x - 
a ni — OTT2- Note that inclusion of the lighcone variable Xi-\ of the quark 
emitter (prior to emission) into p L (fc 7ri |xJ r _ 1 ) was necessary in order to define 
i/xi-i = 1 - a.i/xi-1 in terms of the four momenta fef . 
We are now able to write down the LO+NLO contribution: 

,L+Nf, \ _ 1 -s f d3k 3 f d 3 k 2 f d 3 h 



3 1 ' ' 3! J 2kl J 2k\ J 2k\ 

X $x -x=a 1 +a 2 +a 3 JJ #e*>|ki|>e*0 P3 +N (h , k 2 , k\) , 

i 

^] (P3 (^7r 3 ) ^7T 2 ) ^7Tl) + P3a(^7T3» ^2) ^1 ) + Psbifasi ^""2' ^1 )) > 
7T 

where 



(21) 



p 3 (k 3 ,k 2 ,k 1 ) = p {k 3 \x 2 ) p (k 2 \x 1 ) p (fa\x ) d^^k^^, 

P3ai k 3, k 2,fa) = P L (h\x 2 ) p N (k 2 ,ki\x ) ^kg^maxClkaMkil), (22) 

max(|k 3 |,|k 2 )>|ki|- 



P3bih, k 2 , fa) = p N {k 3 , k 2 \xi) p L (ki\x$) 



With sharp ordering and no interferences. 
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By examining eq. (13) we identify 



p N (k 2 ,k 1 \x') = b 2 {k 2 ,a 2 /x',k 1 ,a 1 /x') 20| k2 |>| kl |, (23) 

where x' of the emitter rescales a's, but transverse momenta kj are un- 
changec 



13 



Two important points about the above equations: (i) For p^ b , where LO 
spectator k 1 ^ proceeds NLO insertion, the ordering | lc2 1 > |ki| is violated, 
i.e. p^ b contributes to two out of six simplices in the |kj| spacd^J This is 
of course general phenomenon and in the case of more such LO spectators 
kj of the NLO insertion may swap with any one of these proceeding LO 
spectators. This shows that the sum is important and instrumental, 
(ii) The validity of the exact NLO DGLAFp^j evolution of eq. Q requires 
that 

D^ +N (t,x)=e- s a 3 (t-t ) 2 



^ ® ® y£(x) + \fi ® t n (x) + ±o>" ® tf(x) 



(24) 



is true. We have checked by analytical and numerical integration that the 
above is almost valid. The only problem is with integration limits for the 
NLO subintegral. Even for simpler p^ a the integration limit of |ki| (evolu- 



tion starting point) is |ki| > e*°, while in eq. (12) this limit is exactly zero. 
One method to overcome this problem and to assure the exact validity of 
eq. Q is to introduce another substantially lower starting point tu of the 
evolution; in the interval (tojiAf) evolution is performed with LO kernels 



and NLO insertions are activated only within (to,t). The error in eq. (24) 
will be reduced to a mere ~ e* M ~*° <C 1 (i.e. power suppressed). 

3.5. Efficient evaluation of the Monte Carlo weight 
The structure of the MC weight for NLO distribution for n partons, 

/ n/2 \ 

r , ht Yl \ Pn (km , ■ ■ ■ , k nn PNRw(k-Kl 1 ■ ■ ■ ) &7T n ) ) 

w _Pn^_ij^ / , s 

Pn Pn\ K iri: ■ ■ • , K n n ) 

includes sum over number of NLO insertions r and sum over set of {w} 
numbering all possible ways of placing r insertions within the n-rung lad- 
der. For many partons it would take prohibitively long CPU time to eval- 
uate the above sums (especially Yl{n}) f° r ever y MC event. However, as 

13 Factor 2#|k,|>|ki| could be omitted - its sole role is to reduce slightly combinatorics. 

14 This is contrary to and p^ a , which contribute to one simplex in i; variables - in 
other words p^ b corrects for the "sharp ordering" of LO in the ] 1<2 1 — ] ki | region. 
Limited, of course, to our subset of diagrams. 
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already pointed out we may cut at r < 2. Moreover, for a given MC point 
(k ni , . . . , k nn ) only one special permutation 7Tfc with the right ordering con- 
tributes in the LO part ^{tt} Pn- I n the NLO part there is a lot of zero 
contributions as well. We exploit this fact and reorder permutations with 
respect to tt^ using group product relation 7r = 7r k 7t' 

n/2 

W n = l + ^L7ir- h ^J2J2J2PNRu( k {* k *'h>---> k (* k *')rJ ( 26 ) 

PnKK**'- r=1 M {7r , } 

and change summation order. Next we restrict summation over the set of 
all permutations {it'} to a much smaller subset {vr 4 "'} of the permutations 
for which p^ RuJ ^ 0: 



w„ 



n/2 

1 + V7\ PNBwi^H^ > ■ ■ • > V^)J ( 27 ) 



How to organize Y1{-k^} sucn that zero contributions are avoided? In the 
general case this is a non-trivial task. Details of the general solution will 
be published elsewhere. Here let us only show it for n = 3. In this case in 



eq. (21) there are two insertions {lv} = {a, b}, see (al-a2) and (bl-b2) in 



Fig. 11 for a only one permutation enters 7r a = (1, 2, 3), while for b only two 



(out of six) permutations ir = (1,2,3), (2, 1,3) do contribute: 

-ii N I N N p N (k2,h\x ) 

W3 = l+W 3a + W 3b , W 3a = A >i M > 

p L {k 2 \x 1 ) p L (h\x ) 

w N = P N {k 3 ,k 2 \xi) g p N (k^ki\x[ 2 ) p L (k 2 \xo) e 

W3b P L Ch\x 2 ) P L {k 2 \xi) h>Hl p L (h\x 2 ) p L (ki\x ) P L {k 2 \xi) h>Hl ' 

(28) 

where x x 2 = xq — a 2 and shorthand notation k; L = k w k, U = t w k was used. 



3.6. Main numerical result 

In the previous section we have sketched the structure of the LO+NLO 
distribution for the MC with exclusive kernels. The above scheme of NLO 
insertions and the methodology of calculating MC weight was worked out 
and tested for up to two NLO insertions for any number of partons in the 



case of the bremsstrahlung diagrams. In Fig. 12 we demonstrate with the 
high precision numerical tests that the above LO+NLO Monte Carlo scheme 
works well in practice. The agreement of two MC results is within statistical 



IFJPAN-IV-09-3 printed on May 9, 2009 



23 



I NLO PG LAP | 




Fig. 12. Comparison of MC results for NLO DGLAP evolutions modelled in the 
exclusive and inclusive ways. Contributions with r = 1, 2 of NLO insertions shown 
separately. Bremsstrahlung diagrams only. The e-term T'^ is omitted. Total 2 • 10 9 
of MC events. 



MC erroi 16 which for r = 1 is ~ 10 3 and for r = 2 



is 



10 4 , relative to 



total LO+NLO. This is the most important result of this report. 

More on what is in the plot in Fig. 12 (i) Both evolutions are imple- 
mented using the same underlying Markovian LO MCj^J (ii) Terms due to 
£ part of 7-traces are omitted. The same exercise, with roughly the same 



precision tag as in Fig. 12 but with e terms included is shown in Fig. 13 



(iii) MC weights are positive, weight distributions are very reasonable, see 
weight distributions in Fig. 14 (iv) Evolution ranges from lOGeV to ITeV; 
LO Prue-evolution, starting from <5(1 — x) at lGeV and ranging to lOGeV, 
provides initial x-distribution for the LO+NLO continuation, (v) As before, 
only Cp part of gluonstrahlung and non-ruining as are used, (vi) NLO vir- 
tual corrections are omitted. However, the methodology of including them 
is at hand, it was already used to include e terms (Fig. 13). 



Obviously, the statistical MC errors are exaggerated in this and the following plots 
by factor 3 or more due to dense binning in log 10 (x). 
17 It can be put easily on top of non-Markovian constrained MC [8]. 
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Fig. 13. Similar test as in Fig. 12 This time e terms are included. 2 • 10 9 events. 
Absolute value is plotted due to the use of the vertical logarithmic scale. 



4. Conclusions and outlook 

The NLO evolution of the parton distribution functions in QCD is a fun- 
damental tool in the lepton-hadron and hadron-hadron collider data anal- 
ysis. It has been formulated for the inclusive (integrated) PDFs with the 
help of inclusive NLO kernels. We report here on the ongoing work in 
which NLO DGLAP evolution is performed for the exclusive (fully uninte- 
grated) multiparton distributions, ePDF's, with the help of the exclusive 
kernels. These kernels are calculated within the two-parton phase space 
using Curci-Furmanski-Petrionzio factorization scheme for bremsstrahlung 
subset of the Feynman diagrams of the non-singlet evolution kernel. The 
multiparton NL O-compatible distribution (with multiple use of the exclu- 
sive NLO kernels) is implemented in the Monte Carlo program simulating 
multi-gluon emission from a single quark emitter. High precision MC results 
show that the new scheme works perfectly well in practice and is fully com- 
patible (equivalent at the inclusive level) with the traditional inclusive NLO 
DGLAP evolution. Once completed, this Monte Carlo module is aimed as a 
building block for the NLO parton shower Monte Carlo, for W/Z production 
at LHC and for ep scattering, and as a starting point for other perturbative 
QCD based Monte Carlo projects. Presented study is limited to Cj, part of 
bremsstrahlung diagrams. In the immediate future the remaining non-single 
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Fig. 14. Weight distributions from LO+NLO Monte Carlo with inclusive kernels 
(upper left) and exclusive kernels (upper right). The x dependence of the weight 
distribution in the latter case is also shown in the lower 2-dim. plot. 

diagrams will be included and the singlet case will be also worked out. 
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